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Abstract 


The changes of relative permeability and capillary pressure as a function of liquid water phase saturation, two key parameters in two-phase 
PEMFC models, are investigated using a capillary network model incorporating an invasion percolation algorithm with trapping. The two- 
dimensional capillary network accounts for capillary dominated drainage and cluster formation. It is shown that relative permeability is constant 
for low saturation, but follows a power law of saturation for high saturations, with an exponent of about 2.4 that is independent of network size or 
heterogeneity. An increase of the network size and reduction in heterogeneity tend to reduce the relative permeability, and relative permeabilities 
of much less then unity are obtained even for saturations as large as 0.8. Capillary pressure on the other hand does not vary with saturation and 
network size, but is influenced by heterogeneity only. This suggests that regardless of the interface shape and size, the capillaries at the interface 
maintain a constant average radius causing the capillary pressure to remain constant. It is finally shown that with appropriate scaling and for a given 
network heterogeneity, the normalized capillary pressure, single-phase permeability and relative permeability can be deduced for other choices of 
porous medium physical scales without requiring a new set of simulations. 


© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


The transport and accumulation of liquid water is a pacing 
item in PEM fuel cells from the view point of design and oper- 
ation, and gives rise to a variety of challenging multiphase flow 
regimes, which are still not well characterized or understood. In 
this paper we focus on the multiphase flow in the porous gas 
diffusion layer of PEM fuel cells. Recent experimental attempts 
to shed light on two-phase flow in PEM fuel cells have been 
made by Tiiber et al. [1], Nam and Kaviany [2], Pekula et al. [3], 
Litster et al. [4], and Bazylak et al. [5] using visualization tech- 
niques and neutron imaging, but there has been little progress in 
quantifying the processes. A variety of modeling strategies have 
been proposed (see review in Litster and Djilali [6]) to predict 
phase saturation in the gas diffusion layers, including the multi- 
fluid model (Berning and Djilali [7]) and the mixture model 
[8,9]. A critical aspect of these models is the need to prescribe 
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the relative permeabilities (k;,;) of each phase and a constitutive 
relation for capillary pressure (pc) as a function of saturation (s). 
Remarkably, multiphase fuel cell models to date have had to rely 
on empirical correlations obtained from soil and sand samples. 
Gas diffusion layers have, as shown in Fig. 1, a radically differ- 
ent structure, and they are also treated to impart hydrophobicity 
in order to promote water transport. The effective relative per- 
meability and capillary pressure curves are expected to be quite 
different [6]. In this paper we propose to address this issue using 
capillary pore network simulations. 

There is a rich literature on multiphase flows in porous media 
due to their relevance in a large number of engineering appli- 
cations ranging from oil recovery, chemical reactors and heat 
exchangers, to drainage and drying of soils. Continuum mul- 
tiphase models are limited in their ability to predict surface 
heat and mass transfer coefficients, distribution of phases (phase 
saturation), and the formation of “dry patches” [10]. Discrete 
pore network models provide an alternative approach to eluci- 
date transfer phenomena and to evaluate transfer parameters. 
In network models, an actual porous medium is represented 
as a network of pores that are connected by throats [11], and 
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Nomenclature 

A system matrix 

A, B, a, p,n power law parameters 

b force vector 

D tri-diagonal matrix 

f fraction of throats invaded at the outlet boundary 
g throat conductance (m? Pa~! s7!) 
kri relative permeability 

K permeability (m7) 

Ksp single-phase permeability (m7) 

l throat length (m) 

L medium length (m) 

nL network size 

N number of random network realizations 
p pressure (Pa) 

Ap/L pressure drop (Pam7!) 

Pe capillary pressure (Pa) 

q flow rate through one throat (m? s7!) 
Q flow rate (m? s7!) 

r radius (m) 

R diagonal matrix 

s phase saturation 

Vv volume (m?) 

Vp volume of pores (m? ) 

W medium width (m) 

Wxôê medium cross-section area (m7) 
Greek letters 

x network heterogeneity, “max//min — 1 
ô medium depth (third dimension) (m) 
0) porosity 

u dynamic viscosity (Pa s) 

o surface tension (Pa m) 

Subscripts 

av average 

bt breakthrough point, f= 1/ny 

H horizontal 

inl inlet 

I fluid phase 

nex next 

out outlet 

p pore 

pre previous 

sp single-phase 

t throat 

tp terminal point, f= 1 

V vertical 


such models have been used to investigate a number of different 
processes, including two-dimensional creeping flow [12], evap- 
oration and drainage [13—17], diffusion and dispersion [18-20], 
and flow in fractures [21]. 

A number of parameters are required for initial network 
definition, which include coordination number, network het- 
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Fig. 1. Scanning electron microscope of a GDL illustrating random fibre struc- 
ture and PTFE webbing. 


erogeneity, pore and throat shape. The coordination number 
(connectivity) is equal to the ratio of connected pores to the site 
pore. For regular networks, the coordination number is equal to 
four and six for 2D and 3D networks, and a value of around 
four has been reported for 3D simulations of stochastic media 
[22]. In direct measurements of porous structures using micro- 
tomography, similar coordination numbers were found; although 
coordination numbers up to fifteen have been reported [22,23]. 
Another parameter defining the network heterogeneity is the size 
distribution of the pores and throats. Random distributions have 
been used with different distribution laws: beta [24], uniform, 
and normal [25]. In addition, the shape and cross-sectional area 
of the pores and throats can be different [26]. Finally, all these 
parameters influence the distribution of the phases in the network 
and interfacial area [27]. 

The multiphase flow parameters (relative permeability, cap- 
illary pressure) depend on the porous medium properties, and 
consequently, network topology [28]. The multiphase parame- 
ters also depend on the type of process, and the predominance of 
gravity, viscous or capillary forces [29]. For very slow processes 
in which capillary forces dominate, the process follows invasion 
percolation [30], for which the problem formulation was pro- 
posed by Wilkinson and Willemsen [31]. The definition of the 
overall potential for any combination of these three forces and 
how the phases distribute within a porous medium is discussed 
in Prat [32] for all other cases. The relative influence of the three 
forces is characterized by the capillary and the Bond numbers 
[13], from which correlation length (average size of immobile 
fluid clusters) can be calculated [33]. The relative permeability 
is constant when the correlation length is much smaller than the 
domain size. Otherwise, when the domain is small compared to 
the correlation length, the relative permeability follows a power 
law of domain size [34]. 

So far the only attempt to apply pore network models to 
fuel cell gas diffusion media is the work of Nam and Kaviany 
[2] who used this method to determine the effective diffusion 
coefficient as a function of porosity and saturation. However, in 
order to compute the required water saturation profile, Nam and 
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Kaviany relied on the prescription of a cubic variation of the 
relative permeability (A;;) with saturation in conjunction with 
the Leverett function for the capillary pressure (pe). The main 
objective of this work is to use network simulations to establish 
the functional dependence of the relative permeability (k,.;) and 
capillary pressure (pc) on the saturation (s) of an invading phase 
(liquid water), particularly in the context of fuel cell gas diffusion 
electrodes. The determination of such constitutive relations is a 
central issue in the development of reliable models for PEM fuel 
cells where two-phase flow conditions are commonly encoun- 
tered [6,7]. In this work, the focus is on the influence of the 
GDL/porous medium domain properties, such as heterogeneity 
and layer thickness, on (;,;) and (pec) in the limit of the capillary 
dominated flow (slow flows). 


2. Capillary network model 


The two main macroscopic properties used to define a porous 
medium are porosity (¢) and permeability (Ksp), and these can 
be interpreted as storage and momentum transfer properties, 
respectively. Capillary network models exploit these in repre- 
senting the medium as a network of pores and throats. The fluids 
are stored in the pores, where there is no resistance to transfer 
between phases, whereas the transfer resistance is associated 
with the throats, defined as volumeless connections between the 
pores. A schematic of such a two-dimensional network with each 
pore connected to four throats is depicted in Fig. 2. In order to 
investigate fluid flow through a specific porous medium using a 
network model, the macroscopic parameters have to be adjusted 
appropriately. Hence, the capillary network is defined such that 
the total volume of pores (Vp) is equal to the void volume of 
the porous medium. Normalizing with the overall volume of the 
porous medium (V), the porosity is obtained: 


ly 
l 


Similarly, the permeability of the network is set by adjusting 
the throats appropriately. In single-phase flow through either a 
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Fig. 2. Porous medium representation as a network consisting of pores and 
throats. 


porous medium or network, Darcy’s law applies. Given a pres- 
sure drop ((Pinl — Pout)/L) across the network, and a flow rate (Q) 
through a network boundary of area (W x ô), the single-phase 
permeability is defined as, 


K : 2 (2) 
oe Pint — Pou Wô 

where (u) is the fluid viscosity. Since the network is a het- 
erogeneous medium, the single-phase permeability (Ksp) is an 
effective permeability that decreases as the medium becomes 
more heterogeneous [35,36]. Hence, the variation of the single- 
phase permeability can be used as an additional constraint to 
bound pore and throat distribution laws. 

The other important issue is how storage and transport quanti- 
ties (porosity and permeability) are altered due to the presence of 
multiphase flow. In multiphase systems, because transfer resis- 
tances and potentials differ, the phase with the lower potential is 
displaced (replaced) by other phase(s) whose potential is higher. 
Thus, the void volume is occupied with two or more phases, 
and in order to quantify phase content, the phase saturation is 
defined. When two or more phases are present, the transport of 
each phase is reduced compared to single-phase flow. Thus, the 
phase permeability (K;) needs to be calculated with respect to 
each phase. To do so, the laws for single-phase flow are assumed 
to be valid in multiphase flow. For the momentum transfer cor- 
responding to each phase (i), the Darcy law relates the phase 
velocity vector (u;) and the phase pressure gradient (Vp;): 

Ki 
u; = —— V pi (3) 


Hi 


where (i= 1, .. ., np) and (mp) is the number of phases. Using Eq. 
(3), the definition of a phase is preserved; this corresponds to a 
part of the domain having different properties than the remaining 
parts. Here, the viscosity (u;i) of the fluid is distinct for each 
phase. Furthermore, given the distribution of each phase, the 
phase parameters for a multiphase system can be calculated. 
The network models can readily be used for such calculations, 
since it is straightforward to obtain the saturation of each phase 
in the overall medium. 

In Eq. (3), the phase permeability (K;) accounts for the 
presence of both porous media and content of other phases (satu- 
ration), whereas the single-phase permeability (Ksp) represents 
the influence of the porous medium. In order to obtain a per- 
meability that accounts for the influence of the saturation only, 
the relative permeability (kr,;) is defined as a ratio of these two 
permeabilities: 


k= Ki a Qi (25) (4) 
l Ksp Api/Li Q 


where (Ap;/L;) is the pressure drop for an individual phase in 
multiphase flow, and is generally different from the single-phase 
pressure drop (Ap/L). Thus, Eq. (4) implies that the pressure 
field can change due to the presence of multiple phases thereby 
altering (kri). 

Since two phases coexist, a condition at the interface(s) 
between the phases has to be included. In general, a pressure dif- 
ference (p; — pj) arises between phases (i) and (j) at the interface, 


B. Markicevic et al. / Journal of Power Sources 171 (2007) 706-717 709 


i.e. the capillary pressure. As a criterion for phase displacement, 
the pressure difference is compared to the threshold capillary 
pressure (pc), which can be calculated from the Young—Laplace 
equation: 

_ 20 


Pc = (5) 
rt 

where (ø) is the surface tension and (r¢) is the radius of capillary. 
For very slow processes (capillary dominated), a simpler dis- 
placement criterion can be used, i.e. the throat with the largest 
radius at the interface will be invaded by invading fluid. The 
invaded fluid can be trapped in the network (clusters) if the 
invading fluid flow paths merge, form internal loops, or hit the 
no-flow network side. Furthermore, the throat invasion is quasi- 
static, meaning that each throat can be invaded at most once, and 
therefore the invading phase pattern is a static configuration as 
well. This process is referred as invasion percolation with trap- 
ping, as defined in Wilkinson and Willemsen [31], and processes 
such as slow drainage follow this pattern. 

Usually, the invasion is stopped once the invading fluid 
appears at the outlet boundary, where only one outlet throat is 
invaded. This is the breakthrough point of the process. However, 
the invasion can be stopped anytime between breakthrough and 
when all throats at the outlet are invaded (the terminal point). The 
fraction of invaded throats at the outlet (f) is therefore (f= 1/nz) 
and (f= 1) for breakthrough and the terminal points, respectively. 
We distinguish two types of invading fluid flow paths: (i) the flow 
paths that connect the network inlet and outlet, (ii) and the flow 
paths that do not reach the network outlet, as they become part 
of invaded fluid clusters. The invading fluid carrying backbone 
comprises all connecting flow paths, and only this part of the 
invading phase contributes to the momentum transfer. The vari- 
ation of the invading fluid flow paths pattern is caused by: (i) 
invasion stopping point, and (ii) network randomness. Using this 
framework, we investigate how multiphase parameters vary with 
invading phase saturation. 


3. Numerical procedure 


We define a regular square capillary network of size (ny, x nz) 
with randomly generated throat radii. The invading fluid enters 
the network at one side (inlet) and the invaded fluid flows out 
of the opposite side of the network (outlet). There is no flow 
through the other two network sides. In a sequence of discrete 
steps, the invading fluid occupies the throat with the lowest 
potential (largest throat radius) at the interface. Thus, the inter- 
face moves and the flow paths of invading fluid expand. The 
invading fluid pattern for a particular network realization is not 
affected significantly by an invasion stopping point (between 
breakthrough and terminal points), as the flow pattern is already 
developed in the breakthrough point. However, the patterns of 
flow paths and clusters can be very different for distinct ran- 
dom realizations that might cause large variations in multiphase 
flow parameters. Two limiting cases are (i) a network with very 
large cluster size, for which flow paths are less branched, and 
(ii) the opposite case of a network with small clusters and 
highly branched flow paths. Fig. 3(a and b) illustrate these 
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Fig. 3. Invasion of a capillary network for cases of (a) low saturation, and (b) 
high saturation of invading phase. The inlet and outlet are located at the top and 
bottom of the network boundary, respectively. 


two limiting cases in which invading and invaded phases are 
shown with thick black and gray lines, respectively. Here, the 
fluid enters the network at the top boundary, and the outlet 
boundary is at the bottom of the network. Corresponding car- 
rying backbones of the case depicted in Fig. 3(b) is shown in 
Fig. 4. 

In order to calculate the single-phase (Ksp) and the invading 
phase permeability (K;) (see Eq. (2) for single-phase permeabil- 
ity), the pressure solutions within the network (for Ksp) and the 
carrying backbone (for K;) are required. This is obtained from 
the simultaneous solution of the material balance equations over 
all pores within the network/backbone. For each pore in the net- 
work/backbone and corresponding throats, the overall material 
balance can be written in the form: 


C 
gee c=4 (6) 
j=l 


where (c) is the coordination number; in single-phase flow this 
represents the number of throats belonging to each pore (equal to 
four for regular square 2D network). In multiphase flow, because 
the backbone is defined as a static configuration with no phase 
transport across an interface, the coordination number for the 
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Fig. 4. Carrying backbone formed for case of high saturation of invading phase. 


pores at the interface is smaller than four since at least one throat 
is occupied by the invaded fluid. 

The material balance for one throat is given by the flow rate 
(qj) through the throat and is calculated from the throat conduc- 
tance (g;) and the pressure difference (Ap) between the starting 
and ending pores of the throat: 


qj = 8j AP (7) 
where the throat is treated as a capillary of radius (r) and length 
(J) occupied by one phase of viscosity (u). The conductance is 
obtained from Poiseuille’s law: 


4 
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In Fig. 5, a pore (i, j) and four connecting throats (i — 1, j), 
(i+ 1, j), (i,j— 1) and (i, j+ 1) are depicted. Both horizontal (H) 
and vertical (V) throats can be previous (pre) or next (nex) to the 
pore (i, j). Since the direction of the flow for each throat is not 
known a priori, all (qj) are defined as flow into (i, j) pore. Using 
this convention and Eq. (6), the material balance and from there 
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Fig. 5. Network node with one pore and connecting throats. 


the pressure balance for pore (i, j) obeys: 


ShreP—10 + Zhex P+10 + Ne P-01 + SrexP+01 
— (8 re + Shox + 8 pre + Zex) Poo = O (9) 


where the subscripts for pressure (p), flow rate (q) and conduc- 
tance (g) are with reference to Fig. 5. An expression of a similar 
form to Eq. (9) can be obtained for a pore at the interface that 
belongs to the carrying backbone, where the throats with invaded 
fluid are omitted as they do not contribute to the balance of the 
invading phase. Applying Eq. (9) to each pore yields a linear 
system of algebraic equations for pressure (p): 


Ap=b (10) 


The vector (b) depends on the pressure boundary conditions, 
and is assembled using Eq. (9) with specified pressures at the 
inlet and outlet of the network. Matrix (A) is a typical finite 
difference matrix and for a network of (ny x ng = n?) pores, 
it has a block-matrix form: 


Rit Dig 0 ... 0 0 0 
Do; R2 Do3 ... 0 0 0 
0 Di R33... 0 0 0 
A=[i i ii: © | aD 
0 0 0 Rn-2,n-2 Dn-2,n-1 0 
0 0 0 Dn-i.n-2 Ry-in-1 Dn-1,n 
0 0 0 tee 0 Dn,n-1 Rn,n 


where matrices (R) and (D) are tri-diagonal and diagonal, 
respectively, with conductance of throats incorporated from Eq. 


(9): 


Pe H H H V V H 
Rii = tridiag(8gpre,i» &pre,i Enex, i Spre,i Snex,i? Enex, 12) 


Dij = diag(g ej (E < J) Stexi @ > J) (13) 


Matrix (A) is a sparse matrix with a scarcity of approx- 
imately 1/ ne. Using sparse matrix solvers for a network of 
ny x ng = 100 x 100 (matrix (A) has 108 elements out of which 
10* are non-zero), the inversion of (A) can be obtained in a few 
seconds (MatLab package). Eq. (11) can be readily extended to 
3D networks and networks with variable coordination number. 

A typical pressure profile for the carrying backbones (Fig. 4) 
of the invading fluid is depicted in Fig. 6. The pressure drop is 
not constant, and the influence of network edges (inlet and out- 
let) is clearly illustrated. As the invading fluid meanders along 
paths with less flow resistance, the pressure changes are different 
in the middle part of the network than at the edges. Averaging 
the pressure along a direction perpendicular to the flow, an aver- 
aged pressure profile along the flow direction is calculated and 
shown in Fig. 7 for four different invasion stopping points: break- 
through (f= 1/n,), terminal (f= 1), and two points (f1) and (f2) 
(/nz <fi <f2< 1), where (f) is a fraction of invaded throats at 
the outlet. Fig. 7 also depicts the pressure drop in the middle 
part of network that is caused by the presence of two-phases (cf. 
Eq. (4) and how relative permeability is calculated). Clearly, the 
pressure drop depends on the stopping point (f). 
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100 


Fig. 6. Calculated pressure profile in the carrying backbone. 


4. Results and discussion 


The simulations are performed using a regular square network 
of (nL X nL) pores with coordination number equal to four. Five 
network sizes were used nz = {20, 40, 60, 80, 100}. A network 
consists of a total of 2m; (nz, + 1) cylindrical throats of dimen- 
sion (r, l), and each throat is occupied by one phase only (either 
invading or invaded fluid). The length of the throat was set to 
(1=2 x 107° m). To account for the porous medium heterogene- 
ity, the throat radius was initially set as arandom variable, with an 
average radius ray =4 x 1074 m. The radii (r) are uniformly dis- 
tributed in range (Tmin, Fmax) X 1074m with (mins Fmax) = {G.5, 
4.5), (2, 6), (1, 7), (0.5, 7.5), (0.2, 7.8)} and a heterogeneity 
parameter defined as (X =Fmax/Tmin — 1); thus x = {0.3, 2, 6, 14, 
38}. For each particular combination (nz, x), N= 100 network 
realizations with randomly generated (r) are obtained, and from 
these we calculate average parameters. Using the above val- 
ues of (ray, J), a single-phase permeability of order 107° m? is 
obtained; whereas for a typical gas diffusion layer (Ksp) is in 
the order of 107!? m?. Scale analysis can be used to adjust the 
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Fig. 7. Averaged (direction perpendicular to the fluid flow) pressure profiles 
along the fluid flow direction for different fraction of invaded throats at the 
outlet (f). The pressure drop in the multiphase region is also shown (dashed 
lines). 


network parameters to achieve the required permeability. In a 
network, the single-phase permeability is proportional to rf l 
and for two networks having the same heterogeneity parameter 
(x), the group Ksp x 1/ rá, should be the same irrespective on 
the values (rav, /). The same scaling applies for the phase perme- 
ability (K;), and therefore, the relative permeability should be a 
function of the heterogeneity parameter (x) only, and should not 
depend on (ray, /). The capillary pressure (pe) is also expected 
to vary proportionally with the inverse of the throat radii at the 
interface (1/r), and the group (Pe X ray) should thus produce 
the same results regardless of the value (ray) as long as (x) 
is kept constant. In order to verify the validity of the scaling 
arguments, the numerical simulations are repeated for the net- 
work nz = 60, but with ray =4 x 1075 m, a value close to those 
observed in GDLs, and /=2 x 1074m and all five x={0.3, 2, 
6, 14, 38}. Using these values of (ray, J), yields (Ksp) values that 
are three orders of magnitude lower (i.e. 10-!? m2), and (Dc) 
values an order of magnitude higher. The relative permeability 
should stay unchanged, and this will be discussed in Section 
4.2. 

The invasion percolation algorithm with trapping is used to 
obtain invading phase distribution into the invaded phase, where 
the invasion is stopped in the range from the breakthrough (one 
throat at the outlet is invaded) to the terminal point (all throats at 
the outlet are invaded). For a known phase distribution, the invad- 
ing phase saturation (s), single-phase (Ksp), invading phase (K;), 
invading phase relative permeability (k,,; = Ki/Ksp), and capillary 
pressure (pc) are calculated. From the pressure solution (see 
Fig. 7), it can be observed that the pressure changes linearly in 
the middle part of the network, where all parameters (s, Ksp, 
Ki, kri, Pc) are calculated. The saturation (s) is calculated based 
on the number of pores occupied by invading phase and overall 
number of pores, whereas permeabilities are obtained from the 
computed pressure and the phase flow rate in the network (Eqs. 
(2) and (4)). Throughout the network, continuity is satisfied and 
the flow rate is constant and equal to the flow rate at the inlet and 
outlet of the network. Finally, the capillary pressure (pc) is found 
by averaging the capillary condition (Eq. (5)) at the backbone 
interface. 

For one network random realization, stopping the invasion 
in the range from breakthrough (f= 1/nz) to the terminal point 
(f= 1), where we recall (f) is the fraction of the invaded throats at 
the outlet, we obtain the invading phase saturation range s = (Sbp, 
Stp). This range changes for distinct random realizations (N). 
Overall, there is an inherently large scatter in the saturations 
obtained from the various realizations (these correspond to the 
local variations in the water content in the GDL). The scatter 
is influenced by both stopping point f=(1/nz, 1) and random 
realization, and is dominated by the latter. Furthermore, the net- 
work size (nz) influences the saturation scatter. For each (nz), 
we calculate the average saturation (Say) for (N = 100) and cor- 
relate this value with (nz) for particular (f), or for all (f). The 
log-log (say—nz) plots are shown in Fig. 8, where (f) is given 
as a parameter. Wilkinson and Willemsen [31] have shown that 
(Say—n) follows a power law: 


Say = Asn“ (14) 
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Fig. 8. Dependence of the average saturation (averaged over As = 1.0) on the 
network size (nz) for different fractions of invaded throats at the outlet (f). 


where the exponent (œ) for invasion percolation with trap- 
ping is (æ =0.18). In our case this corresponds to the stopping 
in the breakthrough point (s= Sbp, f= L/nL), where we found 
(œ =0.183). As we proceed toward (s = stp, f= 1), (œ) decreases 
and approaches (œ =0.10), which is close to the percolation 
exponent (œ =0.12). For (s) averaged for all (f), we found 
(a=0.101). Ultimately, we find that the distribution of (s) does 
not depend on the heterogeneity (x), as spreading of the invading 
fluid is influenced by radius size difference rather than magni- 
tude difference. 


4.1. Single-phase flow 


For single-phase flow, the single-phase permeability (Ksp) is 
computed. Since the throat radius is arandom variable with a par- 
ticular distribution (uniform in this study), the calculated (Ksp) 
is an effective permeability for which it is assumed that Darcy’s 
law applies. The influence of both (nz) and (x) on (Ksp) is investi- 
gated. Providing that the network is stochastically homogeneous 
(sufficiently large), an average of (Ksp) should not depend on 
(nz) (see, e.g. [35,36]), as moments of the throats distribution 
function remain the same. The deviation of (Ksp) is expected to 
decrease as a function of (nz), as the number of throats increase 
with (nz). The results of calculations for (Ksp ~ ny) are sum- 
marized in Fig. 9, where the average (Ksp) and deviation are 
shown with symbols and error bars, respectively. For (nz > 40), 
the average (Ksp) becomes essentially constant. 

As the medium heterogeneity is quantified by the throat dis- 
tribution moments, the effective permeability decreases as the 
porous medium becomes more heterogeneous (see, e.g. [36]). 
According to Katz and Thompson [37], (Ksp) is proportional to 
the characteristic length. The characteristic length is defined with 
respect to some threshold (capillary pressure) in which an infi- 
nite cluster is formed, and this characteristic length decreases 
as medium heterogeneity increases, and thus (Ksp) decreases 
with (x). This implies that (Ksp) is the largest for homoge- 
neous porous media (the corresponding network would consist 
of throats of equal size). It should be noted here that the perme- 
ability of a porous medium is generally at least a two-parameter 
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Fig. 9. Influence of network size (nz) on single-phase permeability. 


function; in the Kozeny equation the two parameters used are 
medium porosity and hydraulic diameter. Katz and Thompson 
[37] showed that the permeability can also be expressed as a 
function of formation factor and a characteristic length with- 
out an explicit dependence on porosity. Similarly, pore network 
simulations appear to yield single-phase permeabilities that are 
independent of porosity, but it should be kept in mind that a 
network is constructed such that its momentum transport resis- 
tance is the same as that of the actual medium. In this work, 
two parameters are again used (fay, l) and these are adjusted 
to achieve the same (Ksp) as the porous medium. However, for 
distinct network realizations, there is a variation in (Ksp). The 
variation of (Ksp) for N= 100 random realizations is shown in 
Fig. 10 for five different (x) and nz = 60. For a network that is 
nearly homogeneous (x ~ 0.3), the variation in (Ksp) is minimal. 
As (x) increases, a larger variation in (Ksp) can be observed, and 
furthermore the average (Ksp) decreases as depicted in the inset 
figure (Ksp ~ x). In Fig. 10, single-phase permeability is plotted 
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Fig. 10. Variations of the single-phase permeability as a function of the het- 
erogeneity parameter for (nz = 60). Inset figure shows (N = 100) realizations for 
(x = 2) and (x = 38). Calculations performed for ray =4 x 1074 and 4 x 1075 m. 
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in terms of the group Ksp x // ae with open and filled symbols 


representing results for ray =4 x 1074m and ray =4 x 1075m, 
respectively. The results demonstrate that scaling by // ri, does 
indeed collapse the (Ksp) results on a single curve as long as 
the heterogeneity parameter is kept constant. The changes in the 
average value and variability of (Ksp) also influence the value and 
variability of the computed relative permeability (k;,;) examined 
next. 


4.2. Relative permeability 


Like saturation (s), the relative permeability of the invading 
phase (k;,;) depends on the network size (nz) and fraction of 
throats invaded at the outlet boundary (f), and there is some (k;;) 
scatter caused by varying stopping points (from breakthrough 
to terminal point) and random network configurations. In the 
context of a gas diffusion layer, the stopping point is particularly 
important due to the presence of the collector plate land area; it 
is only once the liquid water reaches the outlet boundary under 
the channel area that liquid water flow occurs. Regardless of the 
set value of (f), the average of (k;;) decreases with increasing 
(ny). Correlating averaged (,;) and (nz) for a particular (f) we 
find that the power law applies: 


(kni) = Ang’ (15) 


with the exponent (6 = 1.08) for breakthrough point (f= 1/nz) 
and (8=0.914) for terminal point (f=1). Finally, for all (f) 
a value of (8=0.910) is found. These results are shown in 
Fig. 1 1(a). In Eq. (15), we choose (f) and then average the relative 
permeability of invading phase (kr,;) for all random (kr,;) irrespec- 
tive of (s) (this provides a saturation window of As= 1.0). From 
there we deduce the exponent (£). In order to check whether (8) 
depend on saturation (s), we correlate the average (kr) with (nz) 
for a saturation window (As < 1.0). Average (k;,;) versus satura- 
tion (s) is plotted in Fig. 11(b), showing a power law dependence 
of (kri ~ nL) with an exponent (£) that decreases with (s) (in this 
study from 1.05 to 0.719). For low (s), (£) is close to the value 
obtained for (As=1.0), whereas for high (s) there is a large 
decrease in (8). 

Variations in: (i) network size (nz), (ii) stopping point 
(Sbt < S < Stp), and (iii) network configuration (random realiza- 
tions) which correspond to: (i) gas diffusion layer thickness, 
Gi) shoulder/land width, and (iii) local variation in liquid water 
produced, respectively, all generate changes in (kj) and (s). 
These changes are used to investigate the dependence between 
(k,i) and (s), with (k,;) taken as a configurational permeabil- 
ity. In any case (ki) is expected to increase with saturation, 
and for single-phase flow (s=1) it takes the value (k,;= 1). 
Fig. 12 shows how the relative permeability (k;,;) changes with 
saturation (s) for a network size (ny = 100) and a heterogeneity 
parameter (X = max/min — 1 = 38). Each line (kr, ; ~ s) represents 
one random network realization (V=100) for the saturation 
range (Spt < S < Stp); for clarity only some of the (N = 100) real- 
izations are shown. Due to the configuration, different cases arise 
from small to large changes of (s) that are accompanied by large 
or small changes of (k,;) (lines (a) and (b), respectively). Fur- 
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Fig. 11. Dependence of the average relative permeability on the network size 
(nL): (a) curves averaged permeabilities over (As = 1.0) for different fraction of 
invaded throats at the outlet (f); (b) averaged for saturation range (As < 1.0). 


thermore, both (s) and (kr ;) can change significantly (line (c)) 
or exhibit more complex dependence (line (d)). Hence, there is 
a large spread of individual (kr ; ~ s) curves due to the random 
configurations. 
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Fig. 12. Random realization tendencies of (kr, ~ s) curves for stopping points 
range between the breakthrough and the terminal point. 
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Fig. 13. Averaged values of the relative permeability (k;,;) as a function of satu- 
ration: (a) for network random realizations (N = 100) and (nz = 100) and (x = 38) 
in breakthrough point (f= 1/nz), (b) curves (kr; ~ s) for breakthrough, terminal 
point and whole range (1/nz, 1), (c) family of (kr, i ~ s) curves for different (nz). 
Dashed line (in part (c)) represents curve for power of (n=2.4). 


For (ny = 100) and (x =38) at breakthrough point (f= 1/nz), 
the results for (k;.;) and N = 100 random realizations are depicted 
with circles in Fig. 13(a). Invading phase relative permeabil- 
ity (kri) values averaged in saturation increments (here for ten 
increments of As) are depicted with squares, and the standard 
deviation of (k,;) is shown with error bars. For low (s), (k;,i) 


is constant and then increases for higher (s). The procedure 
of averaging (kri) and (s) is repeated for: (i) terminal point 
(f= 1), and (ii) over all range f=(1/nz, 1). Thus, three distinct 
curves (k;,; ~ s) are obtained and shown in Fig. 13(b). All curves 
collapse on essentially the same dependence. Using averaged 
curves (kri ~s) for all (f) (circles and double dotted line in 
Fig. 13(b)) we obtain a family of curves (kr, ~ s) for various 
network sizes (nz) shown in Fig. 13(c). As suggested in the lit- 
erature, there is a power law dependence (k,; ~s”), with the 
power (n) varying from unity to some higher values, depending 
on the capillary number and type of process [38]. For larger net- 
work sizes (nz), (kri < 1) even for high s ~ 0.8 (see Fig. 13(c)). 
Hence, for the non-constant part of a (k;,; ~ s) curve, we correlate 
(k;,i) as a power function of (s): 


kei = Bis” (16a) 


For which (Bz) is found to decrease with (nz), and (n) remains 
constant (n = 2.4). 

On the other hand, by altering the porous medium heterogene- 
ity (x), there is an additional change of (k;.;) However, assuming 
that (;,;) is still a power function of (s), a modified expression 
is adopted: 


kri = Bys" (16b) 


with (n) and (B,) taken to depend on the parameters (nz) and 
(x). The results for constant (nz = 60) and five different (x) are 
presented in Fig. 14. The open symbols in Fig. 14 represent (k;.;) 
results for ray =4 x 1074m, and filled symbols depict results 
for ray =4 x 107° m. The relative permeability shows essentially 
the same dependence (k;,; ~ s) regardless of (ray), consistent with 
the fact that the scaling factor 1/r4, is the same for both (Ksp) 
and (K;) cancels out in the relative permeability (kr, = Ki/Ksp). 
As for (nz), (ki) is constant for low (s) and then increases for 
higher (s). The relative permeability increases with (x) (family 
of curves in Fig. 14). The parameter (B,) increases with (x), and 
an approximately constant value for (n) of around 2.4 is found. 
Here, both powers are equal (Eqs. (16a) and (16b)), regardless 
of (ny) and (x). However, the parameters (Bz) and (B,) and (k;,;) 
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Fig. 14. Relative permeability as a function of saturation for different het- 
erogeneity parameters. Power law fits (n =2.4) are plotted with dashed lines. 
Calculations performed for ray =4 x 1074 and 4 x 1075 m. 
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Fig. 15. Variation of single-phase permeability (Ksp) and invading phase per- 
meability (K;) for (nz =60) with different (x). Errors bars are shown for (Ksp) 
and (K;). (K;) changes more slowly and with larger variation (dash dot circles) 
than (Ksp). 


for low (s) (Fig. 13 and subpart (c, left sides)) vary with (nz) and 
(x). This confirms that it is possible to alter (4;,;) not only with the 
type of process but also with the structure of the porous medium 
(nz, x). However, as the power (n) remains constant, it appears 
that (n) is influenced by the capillary number (process rate) only 
when the correlation length is altered (here the capillary number 
is very low and all changes occur in the capillary regime). 

The most important distinction between the single-phase per- 
meability (Ksp) and the relative permeability (kr,;) is that (kri) 
increases as the heterogeneity parameter (x) increases, i.e. the 
opposite of the behavior of (Ksp).However, the question is how 
does the invading phase permeability (K;) changes with (x)? 
This is addressed in Fig. 15, where (Ksp) and (K;) are corre- 
lated for (nz = 60) at five values of (x). Based on the nature of 
the interaction between the invading fluid and porous medium, 
the external pressure drives the invading fluid into the available 
throats with the largest radius, where the capillary pressure is 
smallest. Hence, as (x) increases, (K;) increases as throats with 
invading fluid become larger, whereas (Ksp) decreases as shown 
in Fig. 10. The results (Ksp, K;) show that (K;) increases with (x) 
more slowly than (Ksp) decreases, but exhibits larger variations 
than (Ksp) (vertical error bars for (K;) are larger than horizontal 
ones for (Ksp)). The error bars reveal a major difference between 
(Ki) and (Ksp), where the variation of (K;) is influenced with 
configuration of the carrying backbone, thus resulting in a large 
variation of (K;), even for small (x). 


4.3. Capillary pressure 


The capillary pressure is calculated at the same part of the cap- 
illary network as the saturation and invading phase permeability. 
At the invaded/invading fluid interface, the pressure difference 
between two neighboring pores, each belonging to two differ- 
ent phases (gas and water), is equal to the capillary pressure, 


Pc = Pnon-wetting — Pwetting- The capillary pressure (pe) is equal to 
the threshold capillary pressure (p. =20/r,) of the throat con- 


necting these two pores. In the threshold capillary pressure, 
the liquid solid contact angle (0) is omitted as in the porous 
medium it may differ from the liquid contact angle at the solid 
surface. If the solid surface contact angle were used, the capillary 
pressure results would shift by a factor cos(@). The additional 
uncertainty in the capillary pressure prediction arises from the 
nature of the gas diffusion layer which can exhibit both wetting 
and non-wetting characteristics. The spreading of the wetting 
phase has been shown to not follow the invasion percolation 
mechanism [33], and networks that account for micro-force bal- 
ance throughout the spreading need to be utilized. The need for 
the micro-force approach arises from the fact that with wetting 
pores present, more than one pore is filled at the same time. 
This violates the basic invasion percolation assumption, hence a 
fully non-wetting medium is modeled in this study. At the inter- 
face, one pore belongs to the cluster of the invaded fluid, and the 
other is filled with the invading fluid. All such pairs of pores are 
identified, and an average capillary pressure is calculated yield- 
ing the capillary pressure in the network. The changes of (pe) 
as a function of both (nz) and (x) are investigated. We observe 
that the capillary pressure behaves similarly to the single-phase 
permeability (Ksp) with respect to both parameters. For vary- 
ing (nL), a constant average (pc) is observed and the deviation 
decreases as (nz) increases (see Fig. 9 for (Ksp)). Similarly, it 
was fond that the stopping point (from breakthrough to terminal 
point that represents the land width) does not influence the cap- 
illary pressure (pe). This constant (pe) is dependent on the type 
of material used for fabricating the gas diffusion layer. 

The variation of (pc) with the invading fluid saturation (s) 
for (nz = 60) and two distinct values of the heterogeneity param- 
eter (x =2 and 38) is shown in Fig. 16 (results for (x =2) are 
given in small inset figure), where pairs (s, pc) are shown with 
symbols. In the same figure, the square symbols with error bars 
show averaged (pc) values and the deviation in saturation incre- 
ments (As). The capillary pressure remains essentially constant 
except for very low (s). This is attributed to the external pres- 
sures at the network inlet and outlet boundaries; these pressures 
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Fig. 16. Discrete random realizations of the capillary pressure for (nz = 60) and 
(x = 38). Inset figure shows the results for (x =2) and the same (nz). 
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Fig. 17. Variation of averaged capillary pressure (symbols) and standard devia- 
tion (bars) as a function of heterogeneity parameter for (nz = 60) for two average 
pore radii values. 


are maintained constant such that the flow is dominated by capil- 
lary forces, and changes of saturations are purely due to random 
network configurations. As the external pressure difference is 
not changed in the present study, the interface shape changes as 
a result of random configurations, but the average throat size at 
the interface remains the same for a set value of (x), and there- 
fore, the capillary pressure remains constant as a function of 
saturation. As (x) increases, the average throat size at the inter- 
face changes, this in turn leads to an increase in (pc) as shown 
in Fig. 17. Two types of symbols (open and full) represent the 
network results for two different networks (ray, /) investigated. 
Scaling the results by (ray), the capillary pressure results pre- 
sented in the form (pe X ray) are independent on value (ray). The 
finding that (pe ~ s) is constant should not be confused with 
(Pe ~ s) curves obtained by varying the external pressure differ- 
ence [39,40], even though almost constant values of (pe) were 
obtained in the middle of the (pe ~ s) curves. When the external 
pressure difference is changed (e.g. increasing the current den- 
sity), the capillary number and, consequently, the interface shape 
are altered, inducing a variation of (pc) with (s). This dependence 
can be expressed using a Leverett type function. Conversely, in 
the present study the capillary number is small and constant, 
thus yielding constant (pe). 


5. Conclusions 


A discrete capillary network model that uses an algorithm 
accounting for invasion percolation with trapping was developed 
and used to investigate the multiphase flow in porous media. 
Numerical simulations were conducted to analyze the behavior 
of two key parameters: relative permeability and capillary pres- 
sure, and their dependence on network size (ny) and network 
heterogeneity (x). The relative permeability decreases as (nz) 
increases, but in contrast with the single-phase permeability, it 
is found to increase with increasing network heterogeneity (x). 
This stems from the fact that the throats with the largest radii are 
preferably invaded, resulting in higher flow rates of the invading 
phase through the network. The relative permeability is found 


to remain constant for low saturation, and to follow a power 
law for higher saturations. The power law exponent does not 
depend on (nz), (x) and stopping point (f) and approximately 
equal to (n =2.4). However, the constant (B) in the relative per- 
meability power law depends on (nz) and (x) suggesting that the 
relative permeability needs to be corrected for the influence of 
the land width (included in parameters (nz)). The simulations 
also suggest that for intermediate saturations (0.2—0.8) the cap- 
illary pressure remains constant for varying (nz) and (f). This is 
attributed to the fact that for a given external pressure, the inter- 
face shape varies, but the average throat size at the interface 
remains the same. Both the relative permeability and capillary 
pressure depend on the heterogeneity (x) which characterizes a 
given porous medium/GDL. Finally, the scaling parameters for 
the relative permeability and the capillary pressure are deter- 
mined and using these scales both (,;) and (pc) are shown to 
reduce to the same dependence irrespective of the network (Fay, 
1), provided that the heterogeneity (x) is kept constant. This 
allows convenient extension of the results to recover parame- 
ters for porous media having different physical scales without 
requiring additional simulations. 
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